Tema 10 - Modelización con `tidymodels`

Técnicas para ‘Big Data’ en Economía - Curso 2023/24
Universidad de Alicante

Pedro Albarrán

Dpto. de Fundamentos del Análisis Económico. Universidad de Alicante

Alberto Pérez Bernabeu

Proceso de modelización

tidymodels es una colección de paquetes para el proceso de modelización (NO implementa modelos) con los principios de tidyverse

install.packages("tidymodels")

library(tidymodels)

Otros paquetes “similares”: mlr3, caret, H2O

Pre-procesado: partición inicial

dplyr transforma los datos para adecuarlos a la modelización, pero tidymodels permite transformaciones específicas para un modelo concreto
initial_split(): particionar los datos en prueba y entrenamiento.

library(mosaicData)
set.seed(9753)
RailTrail_part <- RailTrail %>% initial_split(prop = .8)

censo <- read_csv("data/census.csv") %>% mutate(income = factor(income))
set.seed(7482)
censo_part <- censo %>% initial_split(prop = .8)

Las funciones training() y testing() acceden a cada submuestra

RailTrail_entren <- RailTrail_part %>% training()
RailTrail_prueb  <- RailTrail_part %>% testing()

Pre-procesado: recetas

Las recetas definen las transformaciones a aplicar
recipe() tiene como primer argumento una fórmula (rol de las variables)
Luego se añaden pasos con step_
- step_filter(), step_arrange(), step_rm(), etc.
- step_naomit()
- step_impute_mean(), step_impute_linear(), step_impute_knn()

Se aplican a una variable, todas o un subconjunto con all_outcomes(), all_predictors(), all_numeric(), all_nominal(), contains(), etc.
- step_dummy(all_predictors(), -all_numeric())

Pre-procesado: recetas (cont.)

Se pueden añadir términos polinomiales con step_poly(), interacciones de variables con step_interact(), discretizar variables con step_cut(),etc.
Se centran variables con step_center() o se estandarizan con step_scale()
Se crea un objeto de receta (y luego combinarlo con otras partes del proceso)

receta_lm1 <- training(RailTrail_part) %>%            # Datos: NO crucial
  recipe(volume ~ cloudcover + precip + avgtemp) %>%
  step_poly(avgtemp, degree = 6) %>%                  
  step_center(all_predictors(), -all_nominal()) %>%
  step_scale(all_numeric(), -all_outcomes()) 

receta_logit1 <- training(censo_part) %>%
  recipe(income ~ age + education + race + sex + capital_gain)

También se podría preparar la receta con prep() y aplicar las transformaciones con juice() o bake()

Definir el modelo

tidymodels define un modelo con una interfaz unificada para distintas bibliotecas (con otros nombres de argumentos)

modelo_lm1     <- linear_reg(mode= "regression", penalty = 0) %>%
                    set_engine("lm") 

modelo_glmnet1 <- linear_reg(penalty = 0)  %>% set_mode("regression") %>% 
                        set_engine("glmnet")

modelo_logit1 <- logistic_reg(mode= "classification", penalty = 0) %>% set_engine("glm")

Se podría estimar (entrenar) el modelo con fit()

Flujos de trabajo: `workflow()`

Combina preprocesado y definición del modelo en un único objeto de flujos

lm1_flujo <- workflow() %>%
  add_recipe(receta_lm1) %>%      
  add_model(modelo_lm1)

Un flujo existente se modifica con update_recipe() , update_model(), etc.

Se prepara todo en una única llamada de fit()

lm1_flujo_est <- lm1_flujo %>% fit(data = RailTrail_part %>% training()) 

logit1_flujo_est <-  workflow() %>% add_recipe(receta_logit1) %>%      
                        add_model(modelo_logit1) %>% fit(censo_part %>% training())

Este objeto contiene tanto la receta para transformar los datos como el modelo estimado para mostrar los resultados o predecir

Flujos de trabajo: `workflow()` (cont.)

Se puede extraer la receta para aplicar la transformación a unos datos (p.e., comprobamos que tienen varianza 1 con var())

lm1_flujo_est %>% extract_recipe() %>% bake(RailTrail_part %>% training()) 
lm1_flujo_est %>% extract_recipe() %>% bake(RailTrail_part %>% testing())

También se pueden extraer los resultados de la estimación. Con funciones de broom se convierten a tibbles para trabajar con ellos (p.e., kable() en un documento)

lm1_flujo_est %>% extract_fit_parsnip() %>% tidy()      # resultados de la estimación
lm1_flujo_est %>% extract_fit_parsnip() %>% glance()    # otros detalles de la estimación

broom::augment() calcula predicciones del modelo, residuos, etc.

modelo <- lm1_flujo_est %>% extract_fit_parsnip()
modelo$fit %>% augment()

Predicción

La función predict() (de parsnip) predice valores numéricos, clases, probabilidad de cada categoría, intervalos de confianza, etc.
Devuelve un tibble que podemos añadir a los datos con bind_cols()

 lm1_flujo_est %>% 
  predict(new_data = RailTrail_prueb) %>% 
  bind_cols(RailTrail_prueb %>% select(volume))     # Variable predicha .pred

logit1_flujo_est %>% 
  predict(censo_part %>% testing())                # clase predicha .pred_class

logit1_flujo_est %>% 
  predict(censo_part %>% testing(), type = "prob")    # probabilidades de cada categoría

Notad que se procesan automáticamente los datos donde se predice

Validación del Modelo

Dado un tibble con valores/clases observados (truth) y predichos (estimate), se calcula una métrica (rmse(), rsq(), accuracy(), etc.) o varias (metrics())

lm1_flujo_est %>% predict(RailTrail_prueb) %>% 
  bind_cols(RailTrail_prueb) %>%  
  metrics(truth=volume, estimate= .pred)           #   mae(truth=volume, estimate= .pred)

logit1_flujo_est %>%  predict(censo_part %>% testing()) %>% 
   bind_cols(censo_part %>% testing()) %>%   
  conf_mat(truth=income, estimate= .pred_class)           # clase predicha=
                                                          #  mayor probabilidad predicha

mis_metricas <- metric_set(accuracy, spec, sens)             # métricas específicas
logit1_flujo_est %>%  predict(censo_part %>% testing()) %>% 
  bind_cols(censo_part %>% testing()) %>%
  mis_metricas(truth=income, estimate= .pred_class)

Validación del Modelo (cont.)

Si predecimos probabilidades, se pueden obtener la curva ROC y la AUC

logit1_probs <- logit1_flujo_est %>% 
                    predict(censo_part %>% testing(), type = "prob") %>%
                    bind_cols(censo_part %>% testing()) 

logit1_probs %>% roc_auc(income, `.pred_<=50K`) 

logit1_probs %>% roc_curve(income, `.pred_<=50K`) %>%  autoplot()

Con más de dos clases, se predice la probabilidad de cada clase y la clase predicha es la más frecuente
- La matriz de confusión es similar, pero de dimensiones \(\small k \times k\)
- La accuracy sigue teniendo la misma interpretación
- La ROC-AUC se calcula para cada clase frente a las demás

El proceso de validación cruzada

vfold_cv() crea las particiones en los datos sin procesar

set.seed(9753)
RailTrail_cv <- RailTrail %>% vfold_cv(v=10) # objeto de 10 grupos

Se puede acceder a los datos de entrenamento y prueba de cada uno de ellos con analysis() y assessment(), respectivamente

RailTrail_cv$splits[[1]] %>% analysis() %>% dim()
RailTrail_cv$splits[[1]] %>% assessment() %>% dim()

fit_resamples(), similar a fit(), sobre un flujo de trabajo ya definido y el objeto completo de remuestras de valicación cruzada …

lm1_flujo_cv_est <- lm1_flujo  %>% 
                        fit_resamples(RailTrail_cv)

El proceso de validación cruzada (cont.)

… y el objeto creado contiene los valores de las métricas

lm1_flujo_cv_est %>% collect_metrics()      # promedio sobre 10 iteraciones
lm1_flujo_cv_est$.metrics %>% bind_rows()   # valores en cada iteración

Se pueden cambiar varias opciones del proceso

lm1_flujo_cv_est <- lm1_flujo  %>% 
                        fit_resamples(
                          resamples = RailTrail_cv, 
                          metrics   = metric_set(rmse, mae),
                          control   = control_resamples(save_pred = TRUE)
                        )

Selección de hiperparámetros: tuning

Algunos parámetros, denominados hiperparámetros, no pueden aprenderse directamente durante el entrenamiento del modelo (ej., \(\small \lambda\) en LASSO)

Proceso de ajuste (tuning):
- en una parte de la muestra de entrenamiento, se estiman los parámetros dado un valor del hiper-parámetro
- en la otra parte, se mide el error asociado a ese hiper-parámetro para validar el mejor valor
- se elige el valor con mejor métrica de error en validación

\(\Downarrow\)

\(\hspace{0.1cm}\)

Proceso de tuning

Se pueden identificar los hiperparámetros a ajustar en la receta y/o el modelo

modelo_LASSO <- linear_reg(mode= "regression", 
                           penalty = tune() ) %>%
                    set_engine("glmnet") 

flujo_LASSO_tuning <- workflow() %>%
  add_recipe(receta_lm1) %>% 
  add_model(modelo_LASSO)

flujo_LASSO_tuning %>% parameters()

Se establecen las combinaciones de valores sobre las que se buscará con grid_random(), grid_max_entropy() o grid_regular()

LASSO_grid <- grid_regular(penalty(range = c(0, 15), trans = NULL), levels = 51)
                                   # rango,                 # número de valores

Proceso de tuning (cont.)

Se usa tune_grid() de forma a similar a fit_resamples() en las remuestras por Validación Cruzada de la muestra de entrenamiento

set.seed(9753)
RailTrail_entren_cv <- RailTrail_entren %>% vfold_cv(v=10)

LASSO_tuned <- flujo_LASSO_tuning %>% tune_grid( resamples = RailTrail_entren_cv, 
                                                 metrics   = metric_set(rmse, mae),
                                                 grid      = LASSO_grid               )

Podemos explorar visualmente el ajuste para distintos valores

LASSO_tuned %>% autoplot()

penalty <- LASSO_tuned %>% collect_metrics() %>% filter(.metric == "mae")
penalty %>% ggplot(aes(x=penalty, y=mean)) +  scale_x_log10() +
              geom_line() + geom_point(color="red") + 
              geom_errorbar(aes(ymin=mean-std_err, ymax=mean+std_err), color="gray")

Estimación del modelo completo

NOTA: debemos probar manualmente varios rangos y valores buscando la U
Se puede ver numéricamente el mejor o los cinco mejores candidatos: recordar que, por la variabilidad, varios valores son igualmente aceptables

LASSO_tuned %>% show_best("mae")
mejor_lambda <- LASSO_tuned %>% select_best("rmse")

Actualizamos el flujo de trabajo con un valor, ej., de select_best()

flujo_final <- flujo_LASSO_tuning %>% 
        finalize_workflow(mejor_lambda)  # finalize_workflow(parameters = list(penalty=8.5))

Debemos estimar el modelo en los datos completos de entrenamiento

flujo_final %>%  fit(data = RailTrail_entren) %>%  broom::tidy()

Finalizando y evaluando el modelo

Finalmente, podemos usar last_fit(): ajusta al modelo finalizado en los datos de entrenamiento y lo evalúa en los de prueba.

final_fit <- flujo_final %>% last_fit(RailTrail_part)

final_fit %>% collect_metrics()
final_fit %>% collect_predictions()


final_fit <- flujo_final %>% last_fit(split   = RailTrail_part,
                                      metrics = metric_set(rmse, mae))
final_fit %>% collect_metrics()

Tema 10 - Modelización con tidymodels